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METHOD OF LINES TRANSPOSE: HIGH ORDER L-STABLE 0{N) SCHEMES FOR 
PARABOLIC EQUATIONS USING SUCCESSIVE CONVOLUTION 

MATTHEW F. CAUSLEY*, HANA CHQt, ANDREW J. CHRISTLIEB*, AND DAVID C. SEAL§ 


Abstract. We present a new solver for nonlinear parabolic problems that is L-stable and achieves high order accuracy 
in space and time. The solver is built by first constructing a single-dimensional heat equation solver that uses fast 0{N) 
convolution. This fundamental solver has arbitrary order of accuracy in space, and is based on the use of the Green’s function 
to invert a modified Helmholtz equation. Higher orders of accuracy in time are then constructed through a novel technique 
known as successive convolution (or resolvent expansions). These resolvent expansions facilitate our proofs of stability and 
convergence, and permit us to construct schemes that have provable stiff decay. The multi-dimensional solver is built by 
repeated application of dimensionally split independent fundamental solvers. Finally, we solve nonlinear parabolic problems by 
using the integrating factor method, where we apply the basic scheme to invert linear terms (that look like a heat equation), and 
make use of Hermite-Birkhoff interpolants to integrate the remaining nonlinear terms. Our solver is applied to several linear 
and nonlinear equations including heat, Allen-Cahn, and the Fitzhugh-Nagumo system of equations in one and two dimensions. 

Key words. Method of lines transpose, transverse method of lines, Rothe’s method, parabolic PDFs, implicit methods, 
boundary integral methods, alternating direction implicit methods, ADI schemes, higher order L-stable, multiderivative 


1. Introduction. The prototypical parabolic differential equation is the heat equation. It forms a 
cornerstone of mathematics and physics, and its understanding is essential for defining more complicated 
mathematical models. Fourier introduced this equation as a means to describe transient heat flow. Pick 
quickly recognized its importance to particle and chemical concentrations. As a result, parabolic equations 
are now ubiquitous in describing diffusion processes, which are found in a vast array of physical problems, 
among which are reaction-diffusion models of chemical kinetics [18,19,34], phase field models describing 
morphology and pattern formation in multiphase fluids and solids [3,4,12[, and even the volatility of stocks 
and bonds in mathematical finance [42[ . 

Numerical solutions of (linear and nonlinear) diffusion equations have been the subject of active research 
for many decades. As early as the 1950’s and 60’s, it was recognized that due to the parabolic scaling, method 
of lines discretizations of the heat equation lead to numerically stiff systems of equations, especially for 
explicit time stepping. Larger time steps (on the order of the mesh spacing) can be taken with fully implicit 
solvers, but in practice, full matrix inversions may become difficult and costly, especially when memory is 
extremely limited, as was especially the case for early computers. Thus, alternate dimensionally implicit 
(ADI) splitting methods [11, 13-16,37], that make use of dimensional splitting and tridiagonal solvers, quickly 
gained popularity as part of an effort to reduce the amount of memory required to invert these systems. 

Later on, memory constraints no longer defined the bottleneck for computing, and attention shifted 
toward methods that focused on reducing floating point operations (FLOPs), albeit with additional memory 
constraints. Most notable among these are Krylov methods [22,28], boundary integral methods [20,27], and 
quadrature methods [23,24,30,31,44]. However, with the advent of GPU processors, it appears that we are 
yet again seeing a paradigm shift towards methods that should emphasize small memory footprints, even at 
the expense of incurring a higher operation count. Thus, AD I-like methods, which can efficiently decompose 
larger problems and limit overhead communication, warrant further investigation, and these features are the 
motivating factor for this work. 

In this paper we propose a novel numerical method for obtaining solutions to the linear heat equation, 
and nonlinear reaction-diffusion type equations. As an alternative to classical MOL formulations, we use the 
method of lines transpose (MOL^), which is sometimes referred to as Rothe’s method [38], or the transverse 
method of lines [39] . In this case, the PDF is first discretized in time, and the resulting semi-discrete 
(modified Helmholtz) problem can be solved using a variety of methods. From potential theory [22,27], 
the solution can be constructed by discretizing boundary integrals. However, with dimensional splitting 
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(that is related to the original ADI formulations), the MOL^ can be used to analytically solve simpler, one¬ 
dimensional boundary value problems, and the subsequent solution can be constructed through dimensional 
sweeps, resulting in an 0{N log N) [32,33] or 0{N) [5-7] solver. Furthermore, we extend the method to 
higher orders of accuracy by using a novel idea referred to as successive convolution. This strategy has 
recently been developed in [ 6 ] for the wave equation by the present authors. In the present work, we not 
only extend the method of lines transpose to parabolic problems, but we recognize the resulting expansion as 
a so-called resolvent expansion [ 1 , 21 ], which we leverage to prove stability and convergence of the successive 
convolution series. In addition, we incorporate nonlinear terms with an integrating factor method that relies 
on high order Hermite-Birkhoff interpolants as well as the (linear) resolvent expansions developed in this 
paper. 

The rest of this paper is organized as follows. In Section 2 we derive the basic scheme for the one¬ 
dimensional heat equation, which is L-stable and can achieve high orders of accuracy in space and time. In 
Section 3 we describe how to obtain an arbitrary order discretization in a single dimension with resolvent 
expansions. In Section 4, we describe how this can be extended to multiple dimensions, and in Section 5 we 
present results for linear heat in one and two dimensions. In Section 6 , we describe how our approach can 
handle nonlinear source terms, and in Section 7 we present numerous numerical results including Allen-Cahn 
and the Fitzhugh-Nagumo system of equations. Finally, we draw conclusions and point to future work in 
Section 8 . 

2. First order scheme in one spatial dimension. We begin by forming a semi-discrete solution to 
the ID heat equation using the method of lines transpose (MOL^). Let u = u{x,t) satisfy 


( 2 . 1 ) 


Ut=lUxx, (a;,t) e (a,5) X [0,T], 


with constant diffusion coefficient 7 , and appropriate initial and boundary conditions. The MOL^ amounts 
to employing a finite difference scheme for the time derivative, and collocating the Laplacian term at time 
levels f and t = Following a similar approach from [ 6 ], we introduce a free parameter /? > 0, so 

that the collocation has the form^ 


u"+i - u" 
At 





.n+1 _ 


/3 > 0. 


Next, we introduce the differential operator corresponding to the modified Helmholtz equation, defined by 


( 2 . 2 ) 


L = I- 


d, 

a 


XX 

Y ■ 


p 

VyAf 


After some algebra, we find that the scheme can be written as 


(2.3) 


AK+i - (1 - /32 )w"] = /32u». 


We note that there are at least two reasonable strategies for choosing /3: 

1. Maximize the order of accuracy. For example, if we choose /3^ = 2, then the discretization is 
the trapezoidal rule, which is second order accurate and A-stable. 

2. Enforce stiff decay. For example, if we choose pP = 1, then the discretization is the backward 
Euler scheme, which is first order accurate, L-stable, yet does not maximize the order of accuracy. 

Here and below, we opt for the second strategy, as the stiff decay of numerical solutions of the heat equation 
is of paramount importance. In Section 3.2, we develop this discussion in the context of higher order schemes 
that relies on a careful selection of P as well as repeated applications of a single inverse operator. 

Upon solving equation (2.3) for we find that the equation for the update is 

(2.4) = (1-/?2)u"-^/32/:-1[u"], 


^In [6], there are a total of two time derivatives (and two space), so the right hand side depends on and n" 
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that requires inverting a modified Helmholtz operator. We accomplish this analytically by using Green’s 
functions, from which 

(2.5) = (^ - ^) M := f f e-’^\--y\u{y)dy + . 

The coefficients Ba and Bf, are determined by applying prescribed boundary conditions at x = a, b which we 
describe in Section 2.2. 

Remark 1. Alternatively, had we followed the method of lines (MOL) and first discretized (2.1) in 
space, then the differential operator C would be replaced by an algebraic operator L, and would be inverted 
numerically. 

Remark 2. Although the update (2.4) (with (3 \) is only first order accurate, we describe in Section 

3 how to extend our procedure to arbitrary order in time. 

This MOL^ approach has several advantages. First, the solution is now explicit, but remains uncondi¬ 
tionally stable. Secondly, in recent work [5-7] we show that the convolution integral in Eqn. (2.5) can be 
discretized using a fast 0{N) algorithm, where N is the number of spatial points. We introduce more details 
in Section 2.1, wherein we update the current algorithm to achieve a user-defined accuracy of 0{Ax^) with 
mesh spacing Ax. Finally, since the solution is still continuous in space, we can decouple the spatial and 
temporal errors, and by combining resolvent expansions with dimensional splitting, we extend the method 
to multiple dimensions without recoupling the errors. 

Remark 3. Since dimensional splitting is used, all spatial quantities are computed according to a one¬ 
dimensional convolution integral of the form (2.5), which is performed on a line-by-line basis, following 
so-called "dimensional sweeps". Since the discrete convolution is computed in 0{N) complexity, the full 
solver scales linearly in the number of spatial points (assuming each sweep is performed in parallel). 

A fully discrete scheme is obtained after a spatial discretization of (2.5). The domain [a, b] is partitioned 
into N subdomains [xj-i,Xj], with a = xq < xi < .. .xn = b. The convolution operator is comprised of a 
particular solution, which is defined by the convolution integral 

(2.6) ^M(a:) •= f y e"“l"’"*'lu(y)dy 

and a homogeneous solution 

(2.7) , 

both of which can be constructed in 0{N) operations using fast convolution. We now describe each of these 
in turn, starting with the first. 

2.1. Spatial discretization of the particular solution. The particular solution is first split into 
/[M](a::) = I^{x) + I^{x), where 

I^ix) = I y" e-“("-^)u(y)dy, I^{x) = ^ j' e-“(^-")u(y)dy. 

Each of these parts independently satisfy the first order "initial value problems" 

+ al^{y) = ^u{y), a<y<x, /^(a) = 0, 

=-^u{y), x<y<b, 7^(6) = 0, 

where the prime denotes spatial differentiation. By symmetry, the scheme for 7^ follows from that of 7^, 
which we describe. From the integrating factor method, the integral satisfies the following identity, known 
as exponential recursion 

pXj 

I^{xj) = J^{xj), where = ^ 

2 Jxj.i 
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and 


j — Oi hj^ hj — ^ j ^3 —1 * 

This expression is still exact, and indicates that only the "local" integral needs to be approximated. We 
therefore consider a projection of u{y) onto Pm, the space of polynomials of degree M. A local approximation 

u{xj - zhj) « pj{z), z G [0,1], 

is accurate to 0{h^), and defines a quadrature of the form 

(2.8) = y ^ - hjz)dz ~ y ^ e-''^^pj{z)dz. 

If standard Lagrange interpolation is used, then the polynomials can be factorized as 

r 

(2.9) pj{z) = pjk{z)u,+k = z^Aj^uf, 

k=-l 


where z = [1, z,..., z^]^, and = [uj-g,... ,Uj,..., and Aj is the Vandermonde matrix corre¬ 
sponding to the points xj+k, ioi k = . .r. The values of i and r are such that £ + r = M + 1, and are 

centered about j except near the boundaries, where a one-sided stencil is required. 

Substituting this factorization into (2.8) and integrating against an exponential, we find that 



J^ixj) « j/ : 

= E^ 

where the weights Wj = 

... Wj^r] satisfy 

k=-i 

(2.10) 

T 

= 


and where 

V, H t. 



^3k = H e " z’^dz = 

^ Jo 

2vj 




-E 




p—O 


p\ 


If the weights are pre-computed, then the fast convolution algorithm scales as 0{MN) per time step, and 
achieves a user-defined 0{M) in space. In every example shown in this work, we choose M = 2 or M = 4. 

2.2. Homogeneous solution. The homogeneous solution in (2.7) is used to enforce boundary condi¬ 
tions. We first observe that all dependence on x in the convolution integral, I[x\ := /[u"](a;), in (2.6) is on 
the Green’s function, which is a simple exponential function. Through direct differentiation, we obtain 


( 2 . 11 ) 


h{a)=al{a), h{b) =-al{b). 


Various boundary conditions at a; = a and x = b can be enforced by solving a simple 2x2 system for the 
unknowns Ba and Bf,. For example, for periodic boundary conditions we assume that (at each discrete time 
step, t = f") 


(2.12) u"(a)=u”(6 ), u2{a) = u2{b), Vn G N. 

We next enforce this assumption to hold on the scheme (2.5), 

^“^[^"((a) =/:“^[u”](6) I{a) + Ba + Bi,p = I{b) + BaP + Bb, 

^“^[^"((a) =/:“^[u”](6) a{I{a) - Ba + Bbp) = a{-I{b) - Bap +Bb), 

where p. = and the identities (2.11) are used to find C~^. Solving this linear system yields 


(2.13) 


Ba 


Hb) 

1 - /r’ 


Bb 


Ha) 

1- p' 


Different boundary conditions (e.g. Neumann) follow an analogous procedure that requires solving a simple 
2 x2 linear system for Ba and Bb- 
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3. Higher order schemes from resolvent expansions. In our recent work [6], we apply a successive 
convolution approach to derive high order A-stable solvers for the wave equation. The key idea is to recognize 
the fact that, in view of the modified Helmholtz operator (2.2), the second derivative can be factored as 

(3.1) (“^) = = c-^) ■= CV, 

where 

(3.2) V = I-C-\ C = {I-V)~\ 


Substitution of the second expression into (3.1) determines the second derivative completely in terms of this 
new operator 

(3.3) = 

\ ^ y p=o 


This shows that second order partial derivatives of a sufficiently smooth function u{x) can be approximated 
by truncating a resolvent expansion based on successively applying V to u{x), which is a linear combination 
of successive convolutions (2.5). 

Now, we consider a solution u{x, t) to the heat equation (2.1), that for simplicity we take to be infinitely 
smooth. We perform a Taylor expansion on u{x,t + At), and then use the Cauchy-Kovalevskaya procedure 
[9,40] to exchange temporal and spatial derivatives to yield 


(3.4) 


u{x, t + At) = 


OO 

E 

p—O 


{Atdt 

pi 


-^) = E 

p—O 


{^Atdx 

p\ 


i(x, t) =: 


jAtda: 


i{x, t). 


The term is a spatial pseudo-differential operator, and it compactly expresses the full Taylor 

series. Our goal is to make use of the formula (3.3) to convert the Taylor series into a resolvent expansion. 
This can be performed term-by-term, and requires rearranging a doubly infinite sum. However, if we instead 
work directly with the operator defining the Taylor series, then 


At first glance this expression looks quite unwieldy. However, fortune is on our side, since the generating 
function of the generalized Laguerre polynomials Lp^\z) is 


(3.5) 


p—O 


1 

-p 1-t 

(l-f)^+l 


which bears a striking resemblance to our expansion. Indeed, if we take A = — 1, substitute z = 0^ and 
t = I?, then 


(3.6) = = I + 

p—O p—1 


3.1. Convergence. This expansion has been considered in the context of Cq— semigroups [1,21], where 
(—^^) is replaced with a general closed operator A on a Hilbert space X. In our notation, we restate part 
(ii) of Theorem 4.3 in [1], which is proven therein. 

Theorem 3.1. Let the Cq— semigroup 


= g 4-i)(;32) (j _ 

p=0 ^ ^ ^ 



OO OO 

= ^ {C - If = Y. 

p—O P—O 
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he approximated by 


p 

Tp{p^) = Y,r-^\p)DP. 

p=0 

Then, for u{x) G there exists for each j3‘^ > 0 an integer mg such that for all integers 2 < k < P + 1, 

with P > mo, 

\\T{p)u-Tp{p)u\\ < 

where C{pf‘,k) is a constant that depends only on fT and k. 

Remark 4. The salient point of the theorem is that, in consideration of a (c.f. Eqn. (2.2)), the 
approximation error is of the form CAt^'^^ which matches the form given by a typical Taylor 

method. 

Finally, we truncate the resolvent expansion (3.6) at order p = P. For the heat equation, this defines 
the numerical method as 

p 

(3.7) u{x,t-\- At) = u{x,t) + ^ L^p~^\p^)V^[u]{x,t), 

p^i 

which has a truncation error of the form 

(3.8) T := L^^](/3^)I?'^+^[u](x,<) + 0{At^~'~'^). 

For P = 1, 2,3, these schemes (evaluated at t = t") are 

(3.9) =u"-/32p[u"], 

(3.10) = u” - /32p[u"] - V^[u^], 

(3.11) u"+i = u” - P^V[u^] - V^[u^] - + y) 

We note that for implementation, each operator is applied successively, and is defined by 

p(P+i)[„] := V[VP[u]], V°[u] := u. 

3.2. Stability. There remains one critical issue: the choice of the free parameter [3. In 1974, Nprsett 
studied a similar single-step multiderivative method for the heat equation [35] and he too, had a free param¬ 
eter in his solver. We follow his lead on the Von-Neumann analysis based on his MOL discretization, but in 
this work we optimize (3 to obtain stiff decay, whereas Nprsett chose (3 to maximize the order of accuracy of 
the solver. 

Consider the linear test problem 

^ = Ay, y(0) = 1, A e C, 

whose exact solution y{t) satisfies 

y(t + h) = e^y{t), z = hX G C. 

Application of (3.7) to this test problem results in 

p 

y{t + h) = Y^ {p)bPy{t) = (j){z)y{t) 

p—0 
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where 


D = 




= 1 - (1 - {z/p^)) 


-1 


1 - {zIP^) 

The generalized Laguerre polynomials satisfy many identities, the following of which is the most pertinent: 

X 


(3.12) 


Lpli{x) - L^p\x) = L^^Wx) = 


— (x) 


P +1 J dx 


Here, d'p\x) is the standard Laguerre polynomial Lp{x). Following standard definitions, we say that a 
numerical scheme is A-stable, provided |(/>| < 1 in the left-half of the complex plane z £ C“. Likewise, a 
scheme exhibits stiff decay if (j){z) —>■ 0 as Re{z) —>■ —oo. If an H-stable method also exhibits stiff decay, it 
is L-stable. 

Now, observing that 1 as Re{z) —>■ —oo, we find that 


(3.13) 


lim (l>{z) 

—>■ — 00 • ^ ^ 
p=0 


p=i 


where we have used the first two expressions in (3.12) to introduce a telescoping sum. We are now prepared 
to prove the following. 

Theorem 3.2. Let u{x,t) he an approximate solution to the heat equation (2.1), given by the successive 
convolution scheme (3.7). Then, 

1. If = x\^^ is chosen as the smallest root of L'p_,_j^(x) = {L^pp-^{x)y, then the scheme achieves order 
P -|- 1 , but does not exhibit stiff decay. 

2. If = x^^'^ is chosen as the smallest root of Lp(x) = L^p\x), then the scheme achieves order P, 
and exhibits stiff decay. 

3. Following the first strategy, the schemes are A-stable for P = 1,2,3, whereas the second strategy 
ensures L-stability. For both strategies, A(a)-stability is achieved for P > 3, with large values of 
a « -K12. 

Proof. The proof follows by applying the maximum modulus principle coupled with (3.13). For part 1, 
upon examining the truncation error (3.8), we see that an additional order of accuracy can be gained if we 
choose 

L'pp,{P^) = 0. 


However, Lp{P'^) 0 for this choice, and so stiff decay does not hold. For part 2, we instead enforce 

stiff decay, but then the truncation error is of order P. Finally, part 3 is demonstrated by the maximum 
amplification factors p along the imaginary axis, as shown for both strategies in Figure 1. In particular, we 
observe that \(j){iy)\ < 1 for P = 1,2,3. □ 

Remark 5. In [35], the scheme was chosen to maximize the order of accuracy, implicitly leading to 
eliminating the first term in the truncation error (3.8), which is equivalent to the first strategy. However, in 
this work we follow the second strategy, and choose P^ as the smallest root of Lp{x) to ensure stiff decay. 

For comparison we record the values of /3^ chosen for each order 1 < P < 6 , to those of Nprsett in Table 
1. For all of our solvers, we choose P to be the largest possible value that still yields provable stiff decay. 

4. Resolvent expansions for multiple spatial dimensions. We extend the ID solver to multiple 
spatial dimensions through the use of dimensional splitting. Our key observation is that we can use the 
factorization property of the exponential to perform the series expansion. For instance, in three dimensions, 
we have 

(4.1) el"’- = e-»-(-^)e-'’-(-^)e-»-(-^). 


Now, we first replace each term with the identity (3.6) dimension by dimension, and then truncate the 
expansions which will be in terms of the univariate operators and for 7 = {x,y,zf as defined by 
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(b) Stiff Decay 


Fig. 1: Maximum amplification factors \(j){iy)\ for the first few orders P, with (a) maximal order, or (b) stiff 
decay. When maximizing order, the first 3 schemes exhibit A-stability, whereas ensuring stiff decay leads to 
L-stable schemes. For P > 3, both schemes become A(a)-stable. 


Table 1: Values of /3^ chosen for orders P = 1,2,... 6. The first column are those used in our schemes, and 
uniquely guarantee stiff decay and A(0)-stability. For comparison, we also display the values in Nprsett [35] 
which give optimal order P + 1, at the expense of stiff decay. 



Stiff Decay 

Maximal Order 

p 


Lpm 

/3^ 

Lpm 

1 

1.0000 

0 

2.0000 

-1.0000 

2 

0.5858 

0 

1.2679 

-0.7320 

3 

0.4158 

0 

0.9358 

-0.6304 

4 

0.3225 

0 

0.7433 

-0.5768 

5 

0.2636 

0 

0.6170 

-0.5436 

6 

0.2228 

0 

0.5277 

-0.5211 


(2.2), and (3.7) acting on a function u”(a;, j/, z). This infinite sum with three indices must then be truncated 
to order P, and after a change of indices we find 


(4.2) 


Ep^Y. 


P-1 

P,q,r 


4-i)(/32)p(-i)(/32)L(-i)(/32)i?™; 


y 


in 3D, with the corresponding 2D operator given by 

^P- 1 


(4.3) 


Ep = Y^ 


p,q 


4-i)(/32)p(-i)(/32)ppp«. 


Here we adopt the notation that sums are taken over all non-negative indices that sum to P, and the 
multinomial coefficients are defined such that ( "■ ) = . and ( ” ) = ^. 

\p,q,r/ p\q\r\ P-T- 

Remark 6. The proof of stability for the multi-dimensional algorithm follows direetly from that of the 
one-dimensional case, with the same approach applied to each spatial dimension (i.e. 4>{z) = (j)x{z)4>y{z) for 
the 2D case, and similarly for the 3D case). 
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5. The heat equation. 


5.1. Heat equation in ID. We first illustrate the accuracy of our method for the ID heat equation 
defined in (2.1). We consider initial conditions u(x,0) = sin(a;), for x € [0, 27r] with periodic boundary 
conditions. We integrate up to a final time of T = 4, and set 7 = 0.18^. We use the fast convolution algorithm 
that is fourth order accurate in space (M = 4), and set the spatial grid size to be Ax = « 0.0061. 

This ensures that the dominant error in the solution is temporal. We compute errors by the L°°-norm, and 
compare against the exact solution u{x,T) = e~'^'^uo{x) at T = 4. The result of a temporal refinement 
study for P = 1, 2 and 3 is presented in Table 2. 


Table 2: Refinement study for a ID Heat equation defined in 5.1. 



P = 1 

P = 2 

P = 3 

At 

error 

order 

L°° error 

order 

error 

order 

0.1 

1.8405 X 10-4 

- 

1.6255 X 10-6 

- 

2.4225 X 10-® 

- 

0.05 

9.2121 X lO"*^ 

0.9985 

4.0841 X 10"^ 

1.9928 

3.0620 X 10-y 

2.9839 

0.025 

4.6084 X lO"'^ 

0.9993 

1.0236 X 10"^ 

1.9964 

3.8501 X 10-16 

2.9915 

0.0125 

2.3048 X lO"'^ 

0.9996 

2.5622 X 10-6 

1.9982 

4.8402 X 10-11 

2.9918 

0.00625 

1.1525 X lO"'^ 

0.9998 

6.4097 X 10-6 

1.9990 

6.2021 X 10-12 

2.9642 


5.2. Heat equation in 2D. As a second example, we present results for the 2D heat equation. We 
consider initial conditions u{x,y,0) = sin(x) sin(y), for {x,y) S [0,27r] x [0,27r] with periodic boundary 
conditions. We use a uniform mesh of size Ax = Ay = 2 tt/512 « 0.0123. Likewise, the L°°-error is 
computed by comparing against the exact solution u{x,y,T) — e~‘^~^'^UQ{x,y) at the final time T = 1. In 
Table 3, we present results for a temporal refinement study for orders P = 1, 2, and 3. 


Table 3: Refinement study for a 2D Heat equation defined in 5.2. 



P= 1 

P = 2 

P = 3 

At 

P°°-error 

order 

L°°-error 

order 

^ 00 -error 

order 

0.1 

9.8182 X 10-6 

- 

8.6717 X 10-1 

- 

1.2925 X 10-® 

- 

0.05 

4.9143 X 10-6 

0.9985 

2.1788 X lO-i" 

1.9928 

1.6354 X 10-6 

2.9825 

0.025 

2.4584 X 10-6 

0.9992 

5.4608 X 10-® 

1.9963 

2.0791 X 10-16 

2.9756 

0.0125 

1.2295 X 10-6 

0.9996 

1.3672 X 10-® 

1.9979 

2.9204 X 10-11 

2.8317 


6 . Reaction-diffusion systems. We next extend our method to nonlinear reaction-diffusion systems 
of the form 

(6.1) ut = dv2u + f(u), (x,t) e d X (o,r], 

where u = (ui,U 2 , • • • ,un), with Ui = Ui{x,t), D is a diffusion coefficient matrix, and the reaction term 
F := (/i, / 2 , • • • , /at) is a function of Ui, {i = 1, 2, • • • , N). In the above, C is a bounded domain, 
and we assume appropriate initial values and boundary conditions. We shall view the diffusion as being the 
linear part of the differential operator, and invert this linear part analytically, using successive convolution. 
To derive the scheme, we use operator calculus to first write 

(6.2) (54-DV^)u = F ^ = e-°‘^'F, 
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where e is a pseudodifferential operator. Upon integrating (6.2) over the interval [t,t + At], we arrive 

at the update equation 

/ i+At 

/•At 

(6.3) = / + T)dr, 

Jo 


where we have made use of the abbreviated notation, F(t) := F(u(x, t)). On the left hand side, the diffusion 
terms have been collected by this pseudodifferential operator, and will be approximated using the successive 
convolution techniques developed above. The reaction terms on the right hand side (6.3) are fully nonlinear, 
and we must consider nonlinear stability when choosing a method of discretization. 

We first consider approximating the integral on the right hand side (6.3) with the trapezoidal rule. This 
defines a single-step update equation, which will be second order accurate 


(6.4) 


u{t + At) - 



DAtV^ 


F(t) -h F{t + At) 


We may also obtain a single-step third order scheme, using multiderivative integration [40]. By replacing 
the integrand (6.3) with a third order Hermite-Birkhoff interpolant and performing exact integration of the 
resulting function, we arrive at 


(6.5) u(f + At) - u(t) = 


„DAtV^ 


2 At, 


At 


dF 


— F(t) -k — ( -DAtV^F(t) -f At —(t) 


+ + A^i 


where ^(t) = ^(i) • (DV^u(t) -f F(t)). The Hermite-Birkhoff interpolant that matches the integrand in 
(6.3) at times r = 0, and r = At, as well as its derivative at time t — 0 produces the quadrature rule in 

(6.5) . 

Remark 7. The proposed schemes in (6.4) and (6.5) produce nonlinear equations for u(x,t + At) that 
need to be solved at each time step. Therefore, any implicit solver will necessarily be problem dependent. 

For the problems examined in this work, we make use of simple fixed-point iterative schemes. We 
stabilize our iterative solvers by linearizing F(u) about a background state Fu(u*), which depends on the 
problem under consideration. 

6.1. A discretization of the Laplacian operator. Upon perusing the third order update equation 

(6.5) , we will need to use successive convolution to replace the psuedo-differential operator exp (DAfV^), 
as well as the Laplacian operator V^. This latter point has been detailed in [6], and so we comment briefly 
on it here. Using the one-dimensional expansion (3.3), we observe that the two-dimensional Laplacian is 
similarly given by 


- — = - — _ ^ (j)P I j)P) 

p=i 


and can be truncated at the appropriate accuracy p = P. Here, the subscripts indicate that the convolution 
is only in one spatial direction, and the other variable is held fixed. Thus, Vx is applied along horizontal 
lines for fixed y-values, and likewise for Vy. 

7. Nonlinear numerical results. 

7.1. Allen-Cahn. We examine in greater detail the application of our schem to the Allen-Cahn (AC) 
equation [2], 

(7.1) ut = e'^\'^u + f{u), {x,t) G n X {0,T], 


where the reaction term is f{u) = u — u^, and fl C is a bounded domain, and u satisfies homogeneous 
Neumann boundary conditions. 
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For our fixed point iteration, we linearize / about the stable fixed points u* = ±1, which satisfy 
f'{u*) = 0. For example, the second order scheme from (6.4) becomes 

(7.2) (1 + At) ^ (/"+bfc + 2u"+i’'=) , 

where n indicates the time step as before, and now k is the iteration index. By lagging the nonlinear term 
the fixed point update is made explicit. Likewise, the third order scheme from (6.5) becomes 


(7.3) 



^n+l,fc+l ^ 



2At At 

3 6 


(-e^AtV^/” + Atff) 




2yn+l.fc) _ 


2 yi j 2 

Here, e*^ ^ is again understood by replacing it with a resolvent expansion, which is a truncated series of 

successive convolution operators. 

7.2. Allen-Cahn: One-dimensional test. We demonstrate the accuracy of our proposed schemes 
by simulating a well-known traveling wave solution [8,29], 


(7.4) UAc(a^,^) = ^ - tanh 2^e ^ 0 ) ’ a:G^^ = [0,4], 0<t<T. 

Here, s = ^ = 0.09 is the speed of the traveling wave, and we choose e = 0.03-\/2. Additionally, we choose 
the delay time Tg := 1.5 — sT, so that the exact solution satisfies u^c'(1.5,T) = 0.5. 



Fig. 2: Traveling wave solutions u{x,T) at T = 8 using (7.2) with two different time step sizes, compared 
with the exact profile in (7.4). 


Results for a final time of T = 8 are shown in Figure 2, with two different time steps. The solutions 
agree well with the exact solution. In Table 4, we present the L°°-error in the numerical solution at a final 


Table 4: Refinement study for the ID Allen-Cahn (AC) equation with an exact traveling wave solution 7.2. 



P = 1 

P = 2 

P = 3 

At 

L°° error 

order 

L°° error 

order 

L°° error 

order 

0.025 

2.8216 X 10-^ 

- 

1.3895 X lO"'^ 

- 

2.6060 X 10"*^ 

- 

0.0125 

1.4419 X lO""^ 

0.9686 

3.6115 X lO-t* 

1.9439 

3.9417 X 10"^ 

2.7249 

0.0063 

7.2874 X 10"^ 

0.9845 

9.2164 X 10"^ 

1.9703 

5.5010 X 10-« 

2.8411 

0.0031 

3.6632 X 10"^ 

0.9923 

2.3294 X 10"^ 

1.9842 

7.3122 X 10-9 

2.9113 

0.0016 

1.8365 X 10"^ 

0.9961 

5.8695 X 10-« 

1.9886 

9.5714 X 10-19 

2.9335 


time T = 1, using the exact solution uac{x,T) (7.4). We observe first order accuracy from the Backward 
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Euler method, and the expected orders of accuracy from the second (7.2) and third (7.3) order schemes. To 
ensure that the temporal error is dominant, we have used the fourth order accurate scheme (eq. (2.8) with 
M = 4), with Ax = 2“® to perform spatial integration in the successive convolutions. 

In principle, we can achieve higher orders accuracy in space and time. The latter would require using 
higher order Hermite-Birkhoff interpolation to discretize the reaction term in (6.3). 

7.3. Allen-Cahn: Two-dimensional test. We next solve the Allen-Cahn equation in two spatial 
dimensions. A standard benchmark problem involves the motion of a circular interface [8,29,41], to which 
an exact solution is known in the limiting case e —> 0. The radially symmetric initial conditions are defined 
by 


(7.5) 


i{x, y, 0) = tanh 


' 0.25 - ^(x - 0.5)2 _ 0 5)2 


which has an interfacial circle {u{x,y,0) = 0) centered at (0.5, 0.5), with a radius of Rq = 0.25. This 
interfacial circle is unstable, and will shrink over time, as determined by the mean curvature [2] 


(7.6) 



1 

R' 


Here V is the velocity of the moving interface, and R is the radius of the interfacial circle at time t (i.e., it 
is the radius of the curve defined by u{x,y,t) = 0). By integrating (7.6) with respect to time, we solve for 
the radius as a function of time 


(7.7) 


i?(t) = \/r^- 2eH. 


The location where e is placed in equation (7.1) differs from other references [8,29,41]. Therefore, we point 
out that our time scales have been appropriately rescaled for comparison. 

The moving interface problem was simulated using e = 0.05, At = ^ — = 0.0256, and Aa; = Ay = 

2“® « 0.0039, which are based on the parameters used in [29] . The numerical solution is displayed in Figure 
3, and we observe that the interfacial circle shrinks, as is expected. 



yOOx yOOx yOO 


(a) u{x,y,0) {h)u(x,y,^) (c) u(x,y,T) 

Fig. 3: Time evolution of the initial condition (7.5) of the Allen-Cahn equation, up to time T = 24^^ = 10.24. 

In Figure 4, we plot compare the evolution of the radii obtained by our second order scheme with the 
exact radius (7.7), for two different values of the diffusion parameter e. The radius is measured by taking a 
slice of the solution along y = 0, and then solving for the spatial point where u = 0 using linear interpolation 
between the two closest points that satisfy u(x,0,t) < 0 and u(x,0,t) > 0. Refinement is performed with a 
fixed spatial mesh Ax = Ay = 2“®, and time steps of At = 0.2560,0.1280, and 0.0640. Because the radius 
is derived as an exact solution in the limit (i.e., e —>• 0) [2], we observe that the smaller value of e is indeed 
more accurate. 
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Fig. 4: Radii of the interfacial circle as a function of time (0 < t < ° compared with the reference 
radius R (red line) in (7.7). 


We next perform a refinement study for the Allen-Cahn equations, but this time in two spatial dimen¬ 
sions. To do so, we must incorporate the multivariate successive convolution algorithms in (4.3) and (3.3) 
into the second (7.2) and third (7.3) order schemes. Given that we do not have an exact solution, we compute 
successive errors in an L°°-norm. That is, we compute ||uAt — WAi ||oo for each time step At. Results are as 
expected, and are presented in Table 5. The parameters used are e = 0.05, Aa; = Ay = 2“®, and the final 
computation time is T = 0.5. Again, the quadrature method is fourth order accurate in space, so that the 
dominant source of error is temporal. 


Table 5: Refinement study for the Allen-Cahn equataion in 2D, with homogeneous Neumann boundary 
conditions. 



P = 1 

P = 2 

P = 3 

At 

L°° error 

order 

L°° error 

order 

L°° error 

order 

0.0063 

6.5941 X 10-4 

- 

1.1740 X 10-4 

- 

5.1744 X 10-6 

- 

0.0031 

3.3065 X 10-4 

0.9959 

3.2637 X lO-'^ 

1.8468 

7.8351 X 10-'^ 

2.7234 

0.0016 

1.6563 X 10-4 

0.9973 

8.6726 X 10-6 

1.9120 

1.0811 X 10-^ 

2.8574 

0.0008 

8.2894 X lO"'^ 

0.9987 

2.2389 X 10-6 

1.9537 

1.2961 X 10-« 

3.0602 


7.4. Two-dimensional test: The Fitzhugh-Nagumo system. Finally, we solve a well known re¬ 
action diffusion system that arises in the modeling of neurons, the Fitzhugh-Nagumo (FHN) model [17,25]. 
The FHN system consists of an activator u and an inhibitor v, which are coupled via nonlinear reaction 
diffusion equations 


ut = Du\7'^u+^h{u,v), 

(7.8) d 

Vt = DyV^v + g{u,v), 

where Dy are the diffusion coefficients for u and v, respectively, and 0 < <5 <C 1 is a real parameter. We 
use the classical cubic FHN local dynamics [25] , that are defined as 
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(7.9) 


h{u, v) = Cu{l — u)(u — a) — V, 
g(u, v) = u — dv, 


where C,a and d are dimensionless parameters. The parameters we use are the same as in [10,36]: = 1, 

Dy = 0, a = 0.1, C = 1, d = 0.5, and 6 = 0.005. The diffusion coefficient for the inhibitor is Dy = 0, 
identical to the work found in [26,36,43,45]. 

The second order scheme from (7.2) is applied to each variable u and v separately. This defines the 
numerical scheme as 


(7.10) 


U-+1 = 


= I V 


u 

At 


2S 


+ 


At 


26 


n+1 


where /i” = /i(u",u"') and = g(it"’,u"). We again use a stabilized fixed point iteration to address the 
nonlinear reaction terms. Because {u*,v*) = (0,0) is the only stable excitable fixed point of equation (7.8) 
simply construct the Jacobian of F = {h,g) of (7.9) about this point: 


(7.11) 




1 

* 

s 

1 

s 

1_ 

dih,9)t 

u 


'-Ca 

-l' 


u 


—Cau — V 

1 

* 

1 

_1 


V 


1 

—d 


V 


u — dv 


The resulting second order scheme is 


(7.12) 


CaAt 
1+ 2J 

At ' 


^n+l,k+l 


gAtV= Ln^ 

At 

■1 (/jn+l.fe + yn+l,ky 

At 

- Y 

dAt 


^n+l,fe+l 


At 

yn ^ _gU 

+ 2 

gn+l,k _ yM+l.fc _j_ ^^n+l,fc 


where k is the iteration number, and n is the time level. 

In Figure 5, we present the numerical evolution of the activator u over the domain fl = [—20, 20] x 
[—20, 20], with periodic boundary conditions. We observe similar spiral waves that form in other recent work 
from the literature [10[. 



(a) T = 1 (b) T = 2 (c) T = 4 

Fig. 5: Temporal evolution of the concentration of activator u using (7.12). 


8. Conclusions. In this work, we have introduced a numerical scheme for parabolic problems, which 
achieves high order in space and time for the linear heat equation and some sample nonlinear reaction 
diffusion equations. The scheme is 0{N) for N spatial points, and exhibits stiff decay for any order of 
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accuracy. In the future, we intend to explore parallelization of the algorithm using domain decomposition, 
from which we expect strong performance due to the decoupled spatial factorization that we employ. 

Acknowledgements. We would like to thank the anonymous reviewer and editor for the encourage¬ 
ment to improve the manuscript. This work was supported in part by AFOSR grants FA9550-12-1-0343, 
FA9550-12-1-0455, FA9550-15-1-0282, NSF grant DMS-1418804, New Mexico Consortium grant NMC0155- 
01, and NASA grant NMX15AP39G. 


REFERENCES 

[1] L. Abadias and P. J. Miana. Co-semigroups and resolvent operators approximated by Laguerre expansions. arXiv preprint 

arXiv:1311.7542, 2013. 

[2] S. Allen and J. Cahn. A microscopic theory for antiphase boundary motion and its application to antiphase domain 

coarsening. Acta Metallurgical 1979. 

[3] J. Cahn and J. Hilliard. Spinodal decomposition: A reprise. Acta Metallurgica, 19(2):151 — 161, 1971. 

[4] J. W. Cahn. On spinodal decomposition. Acta Metallurgical 9(9):795 - 801, 1961. 

[5] M. Causley, A. Christlieb, B. Ong, and L. Van Groningen. Method of lines transpose: an implicit solution to the wave 

equation. Math. Comp.., 83(290):2763-2786, 2014. 

[6] M. F. Causley and A. J. Christlieb. Higher order A-stable schemes for the wave equation using a successive convolution 

approach. SIAM J. Numer. Anal., 52(l):220-235, 2014. 

[7| M. F. Causley, A. J. Christlieb, Y. Guclu, and E. Wolf. Method of lines transpose: A fast implicit wave propagator. arXiv 
preprint arXiv:1306.6902, 2013. 

[8] L. Chen and J. Shen. Applications of semi-implicit Fourier-spectral method to phase field equations. Computer Physics 

Communications, 108(2): 147—158, 1998. 

[9] A. J. Christlieb, Y. Giiglii, and D. C. Seal. The Picard integral formulation of weighted essentially non-oscillatory schemes. 

SIAM J. Numer. Anal., 2015. (accepted). 

[10] A. J. Christlieb, Y. Liu, and Z. Xu. High order operator splitting methods based on an integral deferred correction 

framework. Journal of Computational Physics, 294:224-242, 2015. 

[11] J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the 

heat-conduction type. Advances in Computational Mathematics, 6(3):207-226, 1996. 

[12] S. Dai and K. Promislow. Geometric evolution of bilayers under the functionalized Gahn-Hilliard equation. Proceedings 

of the Royal Society A: Mathematical, Physical and Engineering Science, 469(2153), 2013. 

[13] J. Douglas, Jr. On the numerical integration of -h ^ by implicit methods. Journal of the Society for Industrial 

and Applied Mathematics, 3(1):42—65, 1955. 

[14] J. Douglas, Jr. Alternating direction methods for three space variables. Numerische Mathematik, 4(1):41—63, 1962. 

[15] J. Douglas, Jr. and H. Rachford. On the numerical solution of heat conduction problems in two and three space variables. 

Transactions of the American mathematical Society, 82(2):421-439, 1956. 

[16] G. Fairweather and A. Mitchell. A new computational procedure for ADI methods. SIAM Journal on Numerical Analysis, 

4(2), 1967. 

[17] P. C. Fife et al. Mathematical aspects of reacting and diffusing systems. Springer Verlag., 1979. 

[18] R. A. Fisher. The genetical theory of natural selection: a complete variorum edition. Oxford University Press, 1999. 

[19] R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical journal, 1(6):445— 

466, 1961. 

[20] L. Greengard and J. Strain. The fast Gauss transform. SIAM J. Sci. Statist. Comput., 12(l):79-94, 1991. 

[21] V. Grimm and M. Gugat. Approximation of semigroups and related operator functions by resolvent series. SIAM Journal 

on Numerical Analysis, 48(5):1826-1845, 2010. 

[22] J. Jia and J. Huang. Krylov deferred correction accelerated method of lines transpose for parabolic problems. Journal of 

Computational Physics, 227(3):1739-1753, Jan. 2008. 

[23] S. Jiang, L. Greengard, and S. Wang. Efficient sum-of-exponentials approximations for the heat kernel and their applica¬ 

tions. arXiv preprint arXiv:1308.3883, pages 1-23, 2013. 

[24] A. Kassam and L. Trefethen. Fourth-order time-stepping for stiff PDFs. SIAM Journal on Scientific Computing, 

26(4):1214-1233, 2005. 

[25] J. P. Keener and J. Sneyd. Mathematical physiology, volume 1. Springer, 1998. 

[26] V. Krinsky and A. Pumir. Models of defibrillation of cardiac tissue. Chaos: An Interdisciplinary Journal of Nonlinear 

Science, 8(1): 188-203, 1998. 

[27] M. G. A. Kropinski and B. D. Quaife. Fast integral equation methods for the modified Helmholtz equation. J. Comput. 

Phys., 230(2):425-434, 2011. 

[28] J. V. Lambers. Implicitly defined high-order operator splittings for parabolic and hyperbolic variable-coefficient PDE 

using modified moments. Int. J. Pure Appl. Math., 50(2):239-244, 2009. 

[29] H. G. Lee and J.-Y. Lee. A semi-analytical Fourier spectral method for the Allen-Cahn equation. Comput. Math. Appl., 

68(3):174-184, 2014. 

[30] J. Li and L. Greengard. High order accurate methods for the evaluation of layer heat potentials. SIAM Journal on 

Scientific Computing, 31(5):3847—3860, 2009. 

[31] C. Lubich and R. Schneider. Time discretization of parabolic boundary integral equations. Numerische Mathematik, 

455481, 1992. 


15 



[32] M. Lyon and O. Bruno. High-order unconditionally stable FC-AD solvers for general smooth domains II. Elliptic, parabolic 

and hyperbolic PDEs; theoretical considerations. Journal of Computational Physics^ 229(9):3358—3381, 2010. 

[33] M. Lyon and O. P. Bruno. High-order unconditionally stable FC-AD solvers for general smooth domains. H. Elliptic, 

parabolic and hyperbolic PDEs; theoretical considerations. J. Comput. Phys., 229(9):3358-3381, 2010. 

[34] J. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the 

IRE, 50(10):2061-2070, 1962. 

[35] S. P. Nprsett. One-step methods of Hermite type for numerical integration of stiff systems. BIT Numerical Mathematics, 

14, 1974. 

[36] D. Olmos and B. D. Shizgal. Pseudospectral method of solution of the Fitzhugh-Nagumo equation. Math. Comput. 

Simulation, 79(7):2258-2278, 2009. 

[37] D. W. Peaceman and H. H. Rachford, Jr. The numerical solution of parabolic and elliptic differential equations. J. Soc. 

Indust. Appl. Math., 3:28-41, 1955. 

[38] E. Rothe. Zweidimensionale parabolische Randwertaufgaben als Grenzfall eindimensionaler Randwertaufgaben. Math. 

Ann., 102(l):650-670, 1930. 

[39] A. J. Salazar, M. Raydan, and A. Campo. Theoretical analysis of the exponential transversal method of lines for the 

diffusion equation. Numer. Methods Partial Differential Equations, 16(1):30-41, 2000. 

[40] D. C. Seal, Y. Giiglii, and A. J. Christlieb. High-order multiderivative time integrators for hyperbolic conservation laws. 

Journal of Scientific Computing, 60(1):101-140, 2014. 

[41] J. Shen and X. Yang. Numerical approximations of Allen-Cahn and Cahn-Hilliard equations. Discrete Contin. Dun. Sust., 

28(4):1669-1691, 2010. 

[42] S. Shreve. Stochastic Calculus for Einance II: Continuous-Time Models. Number v. 11 in Springer Finance. Springer, 

2004. 

[43] J. Starobin and C. Starmer. Common mechanism links spiral wave meandering and wave-front-obstacle separation. 

Physical Review E, 55(1):1193, 1997. 

[44] J. Tausch. A fast method for solving the heat equation by layer potentials. Journal of Computational Physics, 224(2):956— 

969, June 2007. 

[45] F. Xie, Z. Qu, J. N. Weiss, and A. Garfinkel. Coexistence of multiple spiral waves with independent frequencies in a 

heterogeneous excitable medium. Physical Review E, 63(3):031905, 2001. 


16 



